Exact Hybrid Particle/Population Simulation of Rule-Based Models of Biochemical Systems

Detailed modeling and simulation of biochemical systems is complicated by the problem of combinatorial complexity, an explosion in the number of species and reactions due to myriad protein-protein interactions and post-translational modifications. Rule-based modeling overcomes this problem by representing molecules as structured objects and encoding their interactions as pattern-based rules. This greatly simplifies the process of model specification, avoiding the tedious and error prone task of manually enumerating all species and reactions that can potentially exist in a system. From a simulation perspective, rule-based models can be expanded algorithmically into fully-enumerated reaction networks and simulated using a variety of network-based simulation methods, such as ordinary differential equations or Gillespie's algorithm, provided that the network is not exceedingly large. Alternatively, rule-based models can be simulated directly using particle-based kinetic Monte Carlo methods. This "network-free" approach produces exact stochastic trajectories with a computational cost that is independent of network size. However, memory and run time costs increase with the number of particles, limiting the size of system that can be feasibly simulated. Here, we present a hybrid particle/population simulation method that combines the best attributes of both the network-based and network-free approaches. The method takes as input a rule-based model and a user-specified subset of species to treat as population variables rather than as particles. The model is then transformed by a process of "partial network expansion" into a dynamically equivalent form that can be simulated using a population-adapted network-free simulator. The transformation method has been implemented within the open-source rule-based modeling platform BioNetGen, and resulting hybrid models can be simulated using the particle-based simulator NFsim. Performance tests show that significant memory savings can be achieved using the new approach and a monetary cost analysis provides a practical measure of its utility. © 2014 Hogg et al

HPP performance analyses for various lumping thresholds at cell fraction .

<p>(A) TLBR; (B) Actin; (C) ; (D) EGFR. In all plots, threshold values for different lumping sets are shown on the x-axis. For TLBR and Actin, some thresholds yield the same set of population species as larger thresholds and are thus omitted from the figures. For TLBR, results for thresholds are omitted due to impractically large partial networks in those cases. Results for NFsim (‘NF’) and the hand-picked lumping sets from <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi-1003544-g005" target="_blank">Figs. 5</a>–<a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi-1003544-g008" target="_blank">8</a> (‘HPP’) are shown in all plots for comparison. Error bars show standard error (three samples).</p

Cost of running simulations on the Amazon Elastic Compute Cloud (EC2).

<p>The minimum cost as a function of memory requirement was calculated based on January 2012 pricing (<a href="http://aws.amazon.com/ec2/" target="_blank">http://aws.amazon.com/ec2/</a>) of all <i>Standard</i>, <i>High-CPU</i>, and <i>High-Memory</i> EC2 instances (see Sec. S1 of <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544.s003" target="_blank">Text S1</a> for details of the calculation). Also included are values for NFsim, HPP, and SSA simulations of the EGFR model at cell fraction .</p

Space and time complexities for network-based (SSA) and network-free (NF) stochastic simulation algorithms.

a<p>No dependency graph.</p>b<p>Dependency graph <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Gibson1" target="_blank">[37]</a>, <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Cao1" target="_blank">[38]</a>.</p>c<p>Logarithmic classes (with dependency graph) <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Fricke1" target="_blank">[21]</a>, <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Slepoy1" target="_blank">[39]</a>, <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Schulze1" target="_blank">[40]</a>.</p>d<p>Next-reaction method (with dependency graph) <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Gibson1" target="_blank">[37]</a>.</p>e<p>Direct method (with or without dependency graph) <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Gillespie1" target="_blank">[20]</a>.</p>f<p>Polymerizing systems in gel phase <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Yang1" target="_blank">[23]</a>, <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Monine1" target="_blank">[42]</a> (see <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi-1003544-g005" target="_blank">Fig. 5B</a>).</p>g<p>Direct method-like implementation.</p><p>Scalings are shown with respect to particle number, , and number of reactions, , or rules, . For combinatorially-complex models, . Note that time complexity is given on a “per event” (reaction/rule firing) basis. If a reaction dependency graph <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Gibson1" target="_blank">[37]</a> is used, the space and time complexities of SSA methods with respect to depend on , the maximum number of reactions updated after each reaction firing <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Gibson1" target="_blank">[37]</a>, <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Cao1" target="_blank">[38]</a>. In combinatorially-complex models, often increases with (see <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544.s002" target="_blank">Figure S1</a> of the supporting information). The time complexity of SSA methods with respect to also depends on the method used for selecting the next reaction to fire in the system. Scalings are shown for three different SSA variants that use different selection methods <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Gillespie1" target="_blank">[20]</a>, <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Fricke1" target="_blank">[21]</a>, <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Gibson1" target="_blank">[37]</a>, <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Slepoy1" target="_blank">[39]</a>, <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Schulze1" target="_blank">[40]</a>. Also note that optimized variants of the direct method <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Fricke1" target="_blank">[21]</a>, <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Cao1" target="_blank">[38]</a>, <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-McCollum1" target="_blank">[41]</a> have been shown to outperform methods with lower asymptotic complexity in some cases <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Cao1" target="_blank">[38]</a>. Space and time complexities of the NF algorithm with respect to assume no dependency graph and that the next rule to fire is selected as in Gillespie's direct method <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544-Gillespie1" target="_blank">[20]</a>, although in principle other variants are possible.</p

HPP performance analysis for the EGFR signaling model.

<p>(A) peak memory usage (<i>left</i>: absolute, <i>right</i>: relative to NFsim); (B) CPU run time (<i>left</i>: absolute, <i>right</i>: relative to NFsim); (C) number of reaction events fired during a simulation (); (D) timecourses (means and 5–95% frequency envelopes; ) for activated Sos (<i>top</i>) and nuclear phosphorylated ERK (<i>bottom</i>). The slight deviation from linearity for ‘NF’ in (A) is an artifact of how memory is allocated in NFsim. Due to high computational expense, SSA statistics were not collected in (C) and (D).</p

Partial network expansion (PNE) applied to Rule 11f of Fig. 3.

<p>See Text S4 of the supporting material for the complete, partially-expanded model.</p

Basic workflow of the HPP simulation method.

<p>Given a rule-based model and a user-specified set of population-mapping rules (which define the population species), partial network expansion (PNE) is performed to generate a hybrid version of the original model. The hybrid model is then passed to a population-adapted network-free simulator (e.g., NFsim 1.11), which generates the time-evolution trajectories for all observable quantities specified in the original model.</p

Memory use vs. simulated volume for different simulation methods, including a hypothetical automated HPP (aHPP).

<p>For finite networks, aHPP memory use plateaus once the entire reaction network has been generated. For infinite networks, the scaling at large volumes should fall somewhere between constant and linear (no worse than HPP) depending on the model (see Sec. S2 of <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1003544#pcbi.1003544.s003" target="_blank">Text S1</a> for an analysis).</p

HPP performance analysis for the actin polymerization model.

<p>(A) peak memory usage (<i>left</i>: absolute, <i>right</i>: relative to NFsim); (B) CPU run time (<i>left</i>: absolute, <i>right</i>: relative to NFsim); (C) number of reaction events fired during a simulation (); (D) equilibrium distribution of actin polymer lengths (). The slight deviation from linearity for ‘NF’ in (A) is an artifact of how memory is allocated in NFsim.</p

Simple illustration of ambiguity in the products of reaction rules.

<p>(A) A simple rule encodes the reversible binding of two molecule types, A and B. (B)–(D) If both molecules have multiple binding sites then they may be present within arbitrarily complex complexes. Breaking the bond between A and B thus produces a variety of product species, some of which may correspond to population species and others not. Dashed line represents a bond addition/deletion operation.</p